Solving partial differential equations, using R package ReacTran

Demonstration of the package with plotting in ggplot2 and plotly

Introduction

This notebook is a demonstration of the R package ReacTran. It is a direct implementation of the code and models described in this vignette provided by the package authors, (Soetaert and Meysman). However the plots are implemented using ggplot2, with plotly used for some 3-d examples. The plotting functions included in the package use base R. The purpose of this notebook is to demonstrate the use of modern plotting libraries with the output of the ReacTran package.

Setup

library(ReacTran)
library(tidyverse)
theme_set(theme_bw())
library(patchwork)
library(metR)
library(plotly)

2 A parabolic PDE

As an example of the parabolic type, consider the 1-D diffusion-reaction model, in spherical, cylindrical and cartesian coordinates, defined for \(r\) in [0, 10]: \[ \frac{\partial C}{\partial t} = \frac{\partial}{\partial r}\left(D \cdot \frac{\partial C}{\partial r} \right) - Q \] with \(t\) the time, \(r\) the (radial) distance from the origin, \(Q\), the consumption rate, and with boundary conditions (values at the model edges): \[ \frac{\partial C}{\partial r}_{r=0} = 0 \] \[ C_{r=10} = Cext \]

To solve this model in R, first the 1-D model Grid is defined; it divides 10 cm (L) into 1000 equally-sized boxes (N).

Grid <- setup.grid.1D(N = 1000, L = 10)

Next the properties r and r^2 are defined on this grid:

r <- setup.prop.1D(grid = Grid, func = function(r) r)
r2 <- setup.prop.1D(grid = Grid, func = function(r) r^2)

The model equation includes a transport term, approximated by ReacTran function tran.1D and a consumption term (Q); only the downstream boundary condition, prescribed as a concentration (C.down) needs to be specified, as the zero-gradient at the upstream boundary is the default:

pde1D <- function(t, C, parms, A = 1) {
  tran <- tran.1D(C = C, A = A, D = D, C.down = Cext, dx = Grid)$dC
  
  list(tran - Q)
}

The model parameters are defined:

D <- 1
Q <- 1
Cext <- 20

2.1. Steady-state solution

In a first application, the model is solved to steady-state, which retrieves the condition where the concentrations are invariant, e.g. for the cylindrical coordinate case: \[ 0 = \frac{1}{r^2} \cdot \frac{\partial}{\partial r}\left(r^2 \cdot D \cdot \frac{\partial C}{\partial r} \right) - Q \]

In R , steady-state conditions can be estimated using functions from package rootSolve which implement amongst others a Newton-Raphson algorithm. For 1-dimensional models, steady.1D should be used. The initial “guess” of the steady-state solution (y) is unimportant; here we take simply N random numbers. Argument nspec = 1 informs the solver that only one component is described:

library(rootSolve)
Cartesian <- steady.1D(y = runif(Grid$N), func = pde1D, parms = NULL, nspec = 1, A = 1)
Cylindrical <- steady.1D(y = runif(Grid$N), func = pde1D, parms = NULL, nspec = 1, A = r)
print(system.time(
  Spherical <- steady.1D(y = runif(Grid$N), func = pde1D, parms = NULL, nspec = 1, A = r2)
))

   user  system elapsed 
  0.003   0.000   0.004 

Plot using ggplot2:

data.frame(
  x = Grid$x.mid,
  cartesian = Cartesian$y,
  cylindrical = Cylindrical$y,
  spherical = Spherical$y
) %>% 
  pivot_longer(cols = !contains("x"), names_to = "coord", values_to = "y") %>% 
  ggplot() +
  geom_line(aes(x, y, color = coord, linetype = coord),
            size = 1) +
  scale_x_continuous(breaks = seq(0, 10, 2)) +
  labs(y = "C", color = NULL, linetype = NULL,
       title = "Figure 1: Steady-state solution of the 1-D diffusion-reaction model") +
  theme(panel.grid = element_blank(),
        legend.position = c(1, 0),
        legend.justification = c(1, 0),
        legend.background = element_blank())

2.2. The method of lines

Next the model (for spherical coordinates) is run dynamically for 100 time units using deSolve function ode.1D, and starting with an initially uniform distribution (y = rep(1, Grid$N)):

require(deSolve)
times <- seq(0, 100, by = 1)
system.time(
  out <- ode.1D(y = rep(1, Grid$N), times = times, func = pde1D, parms = NULL, nspec = 1, A = r2)
)

   user  system elapsed 
  0.269   0.000   0.283 

Here, out is a matrix, whose 1st column contains the output times, and the next columns the values of the state variables in the different boxes:

tail(out[, 1:4], n = 3)

       time        1        2        3
 [99,]   98 3.332278 3.332303 3.332366
[100,]   99 3.332383 3.332408 3.332471
[101,]  100 3.332478 3.332503 3.332566

Plot using ggplot2:

# Create a data frame of grid values for joining to the output:
grid_x.mid <- data.frame(
  X = paste0("X", 1:length(Grid$x.mid)),
  grid = Grid$x.mid
)

# Pivot the output to longer format and join with the above data frame to add the Distance values:
out_long <- data.frame(out) %>% 
  pivot_longer(cols = starts_with("X"), names_to = "X", values_to = "density") %>%
  left_join(grid_x.mid, by = "X")

# Plot:
ggplot(out_long) +
  geom_raster(aes(time, grid, fill = density),
             interpolate = TRUE, show.legend = FALSE) +
  geom_contour(aes(time, grid, z = density), color = "black", alpha = 0.5) +
  geom_text_contour(aes(time, grid, z = density), skip = 0, nudge_y = 0.1) +
  scale_x_continuous(breaks = seq(0, 100, 20)) +
  scale_y_continuous(breaks = seq(0, 10, 2)) +
  scale_fill_distiller(palette = "Spectral",
                       breaks = seq(0, 1, 0.2)) +
  labs(x = "Time (days)", y = "Distance (cm)", fill = NULL,
       title = "Figure 2: Dynamic solution of the 1-D diffusion-reaction model") +
  coord_equal(expand = F, ratio = 10)

Add a 3-d plot using plotly:

out_3d <-
  list(
    x = times,
    y = grid_x.mid,
    z = t(out[ , -1])) # need to transpose this to match the x as times and y as r

# Plot:
my_palette <- rev(colorRampPalette(RColorBrewer::brewer.pal(9, "Spectral"))(40)) # ColorBrewer Spectral palette
plot_ly(x = out_3d$x, y = out_3d$y, z = out_3d$z, colors = my_palette) %>%
  add_surface() %>%
  layout(
    scene = list(
      camera=list(
        eye = list(x=0, y=-2, z=1)
        )
    )
  )

3 A hyperbolic PDE

The equation for a wave travelling in one direction (\(x\)) is given by: \[ \frac{\partial ^2 u}{\partial t^2} = c^2 \frac{\partial ^2 u}{\partial x^2} \]

where \(c\) is the propagation speed of the wave, and \(u\) is the variable that changes as the wave passes. This equation is second-order in both \(t\) and \(x\). The wave equation is the prototype of a “hyperbolic” partial differential equation.

For it to be solved in R, the equation is rewritten as two coupled equations, first-order in time:

\[ \frac{\partial u_1}{\partial t} = u_2 \]

\[ \frac{\partial u_2}{\partial t} = c^2 \frac{\partial ^2 u_1}{\partial x^2} \]

We solve the equation with the following initial and boundary conditions: \[ u_1(0, x) = \exp^{-0.05 x^2} \] \[ u_2(0, x) = 0 \] \[ u_{t, - \infty} = 0 \] \[ u_{t, \infty} = 0 \] where the first condition represents a Gaussian pulse.

Here, the grid extends from -100 to 100:

dx <- 0.2
xgrid <- setup.grid.1D(-100, 100, dx.1 = dx)
x <- xgrid$x.mid
N <- xgrid$N

The initial condition, yini and output times are defined next:

uini <- exp(-0.05 * x^2)
vini <- rep(0, N)
yini <- c(uini, vini)
times <- seq (from = 0, to = 50, by = 1)

The wave equation derivative function first extracts, from state variable vector y the two properties u1, u2, both of length N, after which ReacTran function tran.1D performs transport of u1. The squared velocity (\(c^2\)) is taken as 1 (D = 1):

wave <- function (t, y, parms) {
  u1 <- y[1:N]
  u2 <- y[-(1:N)]
  du1 <- u2
  du2 <- tran.1D(C = u1, C.up = 0, C.down = 0, D = 1, dx = xgrid)$dC
  
  return(list(c(du1, du2)))
}

The wave equation can be solved efficiently with a non-stiff solver such as the Runge-Kutta method ode45:

out <- ode.1D(func = wave, y = yini, times = times, parms = NULL, nspec = 2, method = "ode45", dimens = N, names = c("u", "v"))

Plot using ggplot2:

iout <- seq(11, length(times), by = 10)
data.frame(out[c(1, iout), 1:(N+1)]) %>% 
  pivot_longer(cols = contains("X"), names_to = "X", values_to = "u") %>% 
  mutate(x = rep(x, 6)) %>% 
  ggplot() +
  geom_line(aes(x, u, color = factor(time), linetype = factor(time))) +
  scale_x_continuous(breaks = seq(-40, 40, 20)) +
  scale_y_continuous(breaks = seq(0, 1, 0.2)) +
  ggsci::scale_color_d3(name = "time",
                        breaks = c(1, iout) - 1,
                        labels = paste0("t = ", c(1, iout) - 1)) +
  scale_linetype_discrete(name = "time",
                          breaks = c(1, iout) - 1,
                          labels = paste0("t = ", c(1, iout) - 1)) +
  coord_cartesian(xlim = c(-50, 50)) +
  ggtitle("Figure 3: The 1-D wave equation") +
  theme(panel.grid = element_blank())

Plot a 3-d version using plotly:

out_3d <-
  list(
    x = times,
    y = grid_x.mid,
    z = t(out[, 2:(N+1)])) # need to transpose this to match the x as times and y as r

# Plot:
plot_ly(x = out_3d$x, y = out_3d$y, z = out_3d$z, colors = my_palette) %>%
  add_surface()

4 An elliptic PDE

The final example describes again a diffusion-reaction system with production \(p\), consumption \(rC\), and diffusive transport (diffusion coefficients \(Dx\), \(Dy\)) of a substance \(C\) in 2 dimensions \((x, y)\); the boundaries are prescribed as zero-gradient (the default): \[ \frac{\partial C}{\partial t} = \frac{\partial}{\partial x}[D_x \cdot \frac{\partial C}{\partial x}] + \frac{\partial}{\partial y}[D_y \cdot \frac{\partial C}{\partial y}] - rC + p_{xy} \]

The parameter \(p_{xy}\) is the production rate, which is zero everywhere except for 50 randomly positioned spots where it is 1.

Transport is performed by ReacTran function tran.2D; the state variable vector (y) is recast in matrix form (CONC) before it is transported. The first-order consumption rate -r * CONC is added to the rate of change due to transport (Tran$dC). Production, p is added to 50 cells indexed by ii, and which are randomly selected from the grid. The function returns a list, containing the derivatives, as a vector:

pde2D <- function (t, y, parms) {
  CONC <- matrix(nr = n, nc = n, y)
  Tran <- tran.2D(CONC, D.x = Dx, D.y = Dy, dx = dx, dy = dy)
  dCONC <- Tran$dC - r * CONC
  dCONC[ii]<- dCONC[ii] + p
  
  return(list(as.vector(dCONC)))
}

Before running the model, the grid sizes (dx, dx), diffusion coefficients (Dx, Dy), 1st order consumption rate (r) are defined. There are 100 boxes in x- and y direction (n). Furthermore, we assume that the substance is produced in 50 randomly chosen cells (ii) at a constant rate (p):

n <- 100
dy <- dx <- 1
Dy <- Dx <- 2
r <- 0.001
p <- runif(50)
ii <- trunc(cbind(runif(50) * n, runif(50) * n) + 1)

The steady-state is found using rootSolve’s function steady.2D. It takes as arguments amongst others the dimensionality of the problem (dimens) and the length of the work array needed by the solver (lrw = 600000). It takes less than 0.5 seconds to solve this 10000 state variable model.

set.seed(2)
Conc0 <- matrix(nr = n, nc = n, 10)
print(system.time(ST <- steady.2D(y = Conc0, func = pde2D, parms = NULL, dimens = c(n, n), lrw = 6e+05)))

   user  system elapsed 
  0.112   0.000   0.112 

Plot using ggplot:

expand_grid(y = seq(0, 1, length.out = 100), x = seq(0, 1, length.out = 100)) %>% # y first then x to match the output of image
  mutate(z = ST$y) %>% 
  ggplot() +
  geom_raster(aes(x, y, fill = z),
             interpolate = TRUE, show.legend = FALSE) +
  scale_x_continuous(breaks = seq(0, 1, 0.2)) +
  scale_y_continuous(breaks = seq(0, 1, 0.2)) +
  scale_fill_distiller(palette = "Spectral") +
  coord_fixed(expand = F) +
  ggtitle("Figure 5: Steady-state solution of the 2-D diffusion-reaction model") + 
  theme(plot.title = element_text(size = 12))

Plot a 3-d version using plotly:

out_3d <-
  list(
    x = seq(0, 1, length.out = 100),
    y = seq(0, 1, length.out = 100),
    z = matrix(ST$y, nrow = 100, byrow = TRUE)
  )

# Plot:
plot_ly(x = out_3d$x, y = out_3d$y, z = out_3d$z, colors = my_palette) %>%
  add_surface()